function [jac] = fdJac(fhandle,x0,options,exDat0,exDat1,exDat2)
%[jac] = fdJac(fhandle,x0,options,exDat0)
%function to return a finite difference jacobian
%inputs of
%  -function handle for function to finite difference
%  -x0: vector of initial state to perturb about
%  -options: structure that can inclue
%        *.dx0 a scalar or vector of pertubation values for all or each
%        variable
%  -exDat0,1,2: additional inputs to the function

dx0 = 1e-6;

if isfield(options,'dx0')
    dx0 = options.dx0;
end


nExDat = (nargin > 3)*(nargin - 3); % get number of additional data inputs
nVars = length(x0); % number of variables

if length(dx0)~= nVars
    dx0 = ones(1,nVars).*dx0(1);
end

% loop over each variable and perturb system in each direction
for i = 1:nVars
    
    % get variables with ith one perturbed forward
    xp = x0; xp(i) = xp(i) + dx0(i);
    
    % get residual for forward perturbed variables, input differently
    % depending on how many extra data inputs are needed
    switch nExDat
        case 0 
            rp = feval(fhandle,xp);
        case 1
            rp = feval(fhandle,xp,exDat0);
        case 2
            rp = feval(fhandle,xp,exDat0,exDat1);
        case 3
            rp = feval(fhandle,xp,exDat0,exDat1,exDat2);
    end
    
    xm = xp; xm(i) = xm(i) - 2*dx0(i);
    
    switch nExDat
        case 0 
            rm = feval(fhandle,xm);
        case 1
            rm = feval(fhandle,xm,exDat0);
        case 2
            rm = feval(fhandle,xm,exDat0,exDat1);
        case 3
            rm = feval(fhandle,xm,exDat0,exDat1,exDat2);
    end
    
    drdxi = (rp - rm)./(2*dx0(i));
    
    % check to see if r came out as a row or colum vector, flip if it came
    % out as a row to make the jacobian the correct ordering
    [s1,s2] = size(drdxi);
    
    if s2 > s1
        drdxi = drdxi';
    end
    
    
    jac(:,i) = drdxi;
end               